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We employ the macroscopic fluctuation theory to study fluctuations of integrated current in 
one-dimensional lattice gases with a step-like initial density profile. We analytically determine the 
variance of the current fluctuations for a class of diffusive processes with a density-independent 
diffusion coefficient, but otherwise arbitrary. Our calculations rely on a perturbation theory around 
the noiseless hydrodynamic solution. We consider both quenched and annealed types of averaging 
(the initial condition is allowed to fluctuate in the latter situation). The general results for the 
variance are specialized to a few interesting models including the symmetric exclusion process and the 
Kipnis-Marchioro-Presutti model. We also probe large deviations of the current for the symmetric 
exclusion process. This is done by numerically solving the governing equations of the macroscopic 
fluctuation theory using an efficient iteration algorithm. 

PACS numbers: 05.70.Ln, 02.50.-r, 05.40.-a 



I. INTRODUCTION 

Fluctuations around equilibrium states of matter is a 
classical subject of statistical physics. Close to equilib- 
rium, fluctuations of macroscopic observables are fully 
described in terms of the free energy |l| . An important 
recent advance is the elucidation of the behavior of fluctu- 
ations, including large deviations, of macroscopic observ- 
ables in non- equilibrium steady states (NESS) of driven 
lattice gases: simple diffusive transport systems with par- 
ticle conservation pH9|]. The distribution of fluctuations 
in the NESS, as described by the large deviation func- 
tional [lQ], can exhibit qualitatively new features, such 
as non-locality and phase transitions, see Refs. H1G3] 
for reviews. So far, large fluctuations in NESS have only 
been studied in a very few simple lattice gas models. 
These studies, however, have greatly increased the gen- 
eral understanding of fluctuations around NESS. 

Over the last decade a powerful framework, the macro- 
scopic fluctuation theory (MFT) of Bertini, De Sole, 
Gabrielli, Jona-Lasinio, and Landim [135, nas been de- 
veloped for the NESS of diffusive lattice gases driven 
by reservoirs at the boundaries. The MFT is a classical 
Hamiltonian field theory [l3|, [l4{ which describes macro- 
scopic fluctuations in these systems. The MFT formal- 
ism is a further development of the low-noise Freidlin- 
Wentzell theory which in turn is a variant of the 
WKB (after Wentzel, Kramers and Brillouin) approxi- 
mation. A celebrated analog of the MFT for continu- 
ous stochastic systems is the Martin-Siggia-Rose field- 
theoretical formalism [l6[ which has been employed in 
numerous works. Related approaches for lattice gases 
deal, in addition to diffusive transport, with on-site re- 
actions among particles [l7l - [l9| . 

The MFT has been successfully applied to NESS in dif- 
ferent systems [l3|, EH , including those driven not 
from the boundaries. Large fluctuations around NESS 
have also been studied at the microscopic level using ex- 



act [23|-|26| and numerical |27H29j approaches. A perfect 
agreement between the predictions of the MFT and the 
long-time asymptotes of the microscopic calculations has 
been observed whenever the results of the two approaches 
were available. 

With the continuing progress in the studies of NESS, 
a natural next step is to probe fluctuations of macro- 
scopic observables around non- stationary states. Fortu- 
nately, the MFT framework is readily extendable to non- 
stationary settings, such as evolution of a step-like initial 
density profile [30|. There is, however, a major technical 
hurdle which slows down the progress in using the MFT 
for the analysis of both stationary and non-stationary 
problems. Already in the simplest setting of a single 
species of particles, the MFT involves two coupled non- 
linear partial differential equations: the field-theoretical 
Hamilton equations for the density field (a "coordinate" ) 
and a conjugate momentum field. With a few exceptions, 
these equations are not soluble analytically. Still, there 
are several important factors that make the MFT a viable 
alternative to other approaches: 

1. The MFT is stripped off unnecessary details of 
microscopic interactions, so it directly probes the 
large-scale, long-time asymptotic regime that is of 
most interest. 

2. The MFT provides the "optimal path": the density 
profile history which gives a dominant contribution 
to the probability of observing, say, a given current. 

3. The Hamilton equations underlying the MFT can 
be solved numerically with an iteration algorithm 
|3lj . Alternatively, a numerical minimization of 
the mechanical action, intrinsic in the MFT, can 
be performed [321 ] . These numerical algorithms are 
much more computationally efficient than micro- 
scopic stochastic simulations. 
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4. As we show here, a perturbative analytic solution is 
possible which probes, for a whole class of models, 
small fluctuations. 

The main objective of this work is to demonstrate 
these advantages. We will investigate, within the MFT 
formalism, the noisy evolution of a step-like initial den- 
sity profile in a class of lattice gas models in one dimen- 
sion, where the transport is symmetric and diffusion-like. 
Two non-trivial examples are the simple symmetric ex- 
clusion process (SSEP) which has been extensively stud- 
ied (see e.g. [2H9| and references therein) and the Kipnis- 
Marchioro-Presutti (KMP) model [H-[35j]. (In the SSEP 
each particle can hop to a neighboring site at rate 1 if 
that site is unoccupied by another particle. If it is occu- 
pied, the move is forbidden. The KMP model is a one- 
dimensional chain of mechanically uncoupled harmonic 
oscillators which randomly redistribute energy among 
neighbors.) The step-like initial condition for the par- 
ticle density, 

= o) = ('-' x< °> (1) 

[p+, x > 0, 

provides a good "litmus test" for theory of fluctuations 
in non-stationary systems. As in Refs. [13, H3], we will 
be interested in the statistics of integrated current - the 
total number of particles or the total energy — passing 
into the half-line x > during a given time T. The 
precise mathematical formulation of the problem in the 
framework of the MFT was given a few years ago [13] , but 
the problem has defied solution except for the completely 
integrable case of non-interacting random walkers. Our 
strategy here will be to solve the problem perturbatively 
around the noiseless hydrodynamic solution, thus prob- 
ing typical, small fluctuations of the current. In addition, 
we will show how large current fluctuations can be effi- 
ciently simulated numerically. 

We will consider a class of models whose hydrodynamic 
description is provided by a diffusion equation 

d t p = d x [D{p) d xP ] (2) 

with the diffusion coefficient D(p). Having solved this 
equation with the initial condition ([ij, one can compute 

/>oo 

(J(T)) - / dx[p(x,T)-p(x,0)}, (3) 
Jo 

the average integrated current for the underlying micro- 
scopic model. At the level of MFT, the class of micro- 
scopic models that we consider here is fully characterized, 
in addition to the diffusion coefficient D(p), by the func- 
tion cr(p) which describes equilibrium fluctuations 0. 
Our formalism can handle, in a simple way, systems with 
a density-independent diffusion coefficient (which we set 
to D = 1 without loss of generality) but an arbitrary 
a(p). 

The total current J = J(T) into the right half-line is 
a random quantity. The average total current grows as 



VT in the long time limit, T> 1. The variance of the 
total current, (J 2 ) c = (J 2 ) — (J) 2 , also exhibits a diffusive 
growth with time: 

(J 2 ) c = V(p.,p+,a)Vf. (4) 

The quantity V depends on the densities p± and, through 
(7 = <j(p), on the model. Intriguingly, one has to be care- 
ful in defining the averaging procedure [3fJ, |37[ . In the 
quenched setting the initial condition (TTJ) is determinis- 
tic. In the annealed setting one allows equilibrium fluc- 
tuations in the initial condition (JT]). More precisely, the 
initial density profile in the left (correspondingly, right) 
part of the system is chosen from the equilibrium prob- 
ability distribution corresponding to density p- (corre- 
spondingly, p+). As a result, the most probable initial 
density profile, see below, is different from a step func- 
tion. 

The main analytical results of this work are explicit 
expressions for V for diffusive processes with D = 1 and 
arbitrary u{p). In the quenched setting we obtain 

Vquenchcd = / ~. — 7 / dx (J [p(x, I - <)] e^ 2 * , (5) 
JO 47rf J -co 

where p(x, t) is the solution of the classical diffusion equa- 
tion with D = I and initial condition (J]), see Eq. (|36|) 
below. In the annealed setting we obtain 

\/2 — 1 

Knnealcd = ^quenched + /— W(P-) + <?(P+)] • (6) 

2y 2n 

Since o~{p) is intrinsically positive (for p > 0), the second 
term on the right-hand side of Eq. ([6]) is positive. Hence 
Knneaiod > Vq U enched, as expected on physical grounds. 

The rest of this paper is organized as follows. SectionlTTl 
includes important preliminaries that will be used in the 
subsequent sections: We briefly discuss the moment gen- 
erating function of the current and its long-time behavior, 
formally introduce the functions D(p) and <r(p), and out- 
line the MFT formulation, due to Derrida and Gerschen- 
feld [13] , of the problem of statistics of the current for the 
step-like initial condition ([!]). In Sections IIIII and ITVl we 
develop a perturbation theory around the noiseless hy- 
drodynamic solution and determine the variance of cur- 
rent, alongside with the optimal paths, in the quenched 
and annealed settings. Particular examples of these re- 
sults for the symmetric state, p- = p+, for the SSEP and 
KMP models, and for the non-interacting random walk- 
ers, are presented in Section [V] Section PVTl is devoted to 
a numerical calculation, within the MFT, of the optimal 
path conditioned on observing a large deviation of the 
current. Concluding remarks appear in Sec. lVIll Finally, 
in one of the Appendices we present an alternative way 
of calculating the variance in the quenched setting: by 
employing fluctuating hydrodynamics. 
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II. PRELIMINARIES AND GOVERNING 
EQUATIONS 

A. Moment generating function 

A complete description of the current fluctuations is 
provided by the probability distribution P(J,T). Often 
it is more convenient to deal with the moment generating 
function 



^e AJ P(J,T) = 1 



(7) 



,7>0 



that encapsulates P(J,T) and provides the moments of 
the distribution. Alternatively, by taking the logarithm 
of IJ7J, one can rewrite this expression as 



ln(e AJ ) 



A n 



(J") 



(8) 



As an example, Fig. Q] shows two plots for Mannaaicd/v 77 ?" 
versus A: the exact long-time result from Eq. (|10j) and the 
two-cumulant approximation p, — C\\ + (1/2) C2A 2 , for 
p_ = 0.6 and p+ = 0.2. As one can see, a discrepancy 
appears only at |A| ~ 5. Correspondingly, deviations 
of the probability distribution P(J,T) from gaussianity 
only occur in far tails of the distribution. 




which defines the cumulants of the distribution. The first 
two cumulants are (J) c = (J) and (J 2 ) c = (J 2 ) — (J) 2 - 

In diffusive systems with the step-like initial condition 
([T]), the moment generating function exhibits the follow- 
ing long-time behavior: 



3 Vt »(\, P - iP+ ) 



(9) 



For non-interacting random walkers (RWs), the function 
/i(A, p-,p+) was calculated [30] , both in the quenched 
and in the annealed settings, from the MFT formalism. 
For the SSEP it was also calculated (36[, from the micro- 
scopic model, in the annealed setting: 



SSEP _ j_ 

"annealed _ 
7T 



dk In ( 1 + Ae 



-If 



(10) 



FIG. 1: (Color online). Plotted versus A are the exact long 
time result for M^nne^icd/v 77 ? from Eq. (Tl0)l (the solid curve) 
and the two-cumulant approximation CiA + (1/2) C2A 2 with 
Ci and C2 from Eqs. (|12|) and (|13[) (the dashed curve) for 
p_ = 0.6 and p+ = 0.2. The two circles are numerical results 
obtained with iteration algorithm described in Section IVll 



Derrida and Gerschenfeld [30] also found the function 
Mannoaicd(A,p-, / o+) for the KMP model. They showed, 
within the MFT formalism, that it is related to //annealed 
for the SSEP, so it can be established without a new 
calculation: 



KMP 
"annealed 



1 

2^ 



dk In 1 



(l + Ae-* 2 ) 



(14) 



where 
A = P-(e A 



where 



l) + P+(^ A -l)+P-P+(e A -l)(e~ 



1 



In the special case of p- = 1 and p+ — 0, the initial 
condition in the SSEP cannot fluctuate. As a result, 
/i(A,l,0) is the same for both annealed and quenched 
settings. 

Expanding the integrand in Eq. (|10p in powers of A, 
we can extract the cumulants for the SSEP. They have a 
universal long-time behavior 



(JP) C = C>_, P+ )VT 



with 



^C 2 



P- 
P- 
1 



P+, 
P+ - 



P 2 - 



3*) I*- 



P'-\ 



(11) 



(12) 



(13) 



etc. Although obtained via an expansion in small A, the 
cumulants C\ and C2 provide a surprisingly good approx- 
imation of p for |A| comparable to or even greater than 1. 



A = 2A(p+ -/>_)- 4A 2 p_P+. 



(15) 



Expanding the integrand in Eq. (1141) in powers of A, 
yields Eq. (ITTj) . with the same C\ as for the SSEP, and 



Notably, the average current, 

P- - P-i 



(J) 



T 



(16) 



(17) 



is the same for the annealed and quenched averages, and 
for any model with D(p) = 1, including non- interacting 
random walkers, the SSEP and the KMP model. The 
variance is already model-dependent, and it also depends 
on the type of averaging. The above expressions for C 2 
for the SSEP and KMP models refer to the annealed 
case. To our knowledge, in the quenched case even the 
variances are unknown; they will be in the focus of this 
work. 
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Model 


D(r) 


cr(r) 


F(r) 


RW 


1 


2r 


r In r — r 


SSEP 


1 


2r(l - r) 


r lnr + (1 — r) ln(l — r) 


KMP 


1 


4r 2 


-(1/2) lnr 



TABLE I: Functions D(r),a(r),F(r) for non-interacting ran- 
dom walkers and for two interacting particle systems, the 
SSEP and the KMP. 



B. 



D, a and F 



Here is a brief recap of the formal definitions of the 
quantities D(p) and cr(p), and of their relation to the 
free energy density in equilibrium, F(p) Q. The func- 
tions D(p) and u{p) characterize the flux and its variance, 
respectively, in a simple stationary setting. Consider a 
one-dimensional system of a finite (but very large) length 
L which is in contact with reservoirs of particles (or en- 
ergy) with density p- on the left and p+ on the right. 
When these densities are close to each other, p± ~ r, 
with 



\P+ ~P-\<r, 



(18) 



the average flux per unit time is proportional to D(r) 
and the density difference: 



r (J) D(r) 

hm — = —j—{p- -P+). 



(19) 



In its turn, er(r) can be extracted from the growth law 
of the variance of the flux evaluated at the equilibrium 
state p± = r: 



(J 2 



lim 

t— >oo t 



cr(r) 
L ' 



(20) 



Therefore, the quantities D(r) and <j(r) characterize 
small deviations from equilibrium. The equilibrium ori- 
gin of D(r) and cr(r) is additionally emphasized by the 
equation 



d 2 F 
dr 2 



2D(r) 
cr(r) 



(21) 



relating D(r) and cr(r) to the equilibrium free energy 
density F(r). Equation ([21} follows 0, [Tl| from the 
fluctuation-dissipation theorem. It also appears natu- 
rally in the MFT formalism, see Appendix A. 

Table I lists the functions D(r), a(r) and F(r) for three 
specific models: the RW, the SSEP and the KMP. 



C. MFT formalism 

The MFT formalism QJ, Q HI describes large devia- 
tions of macroscopic quantities in diffusive lattice gases. 
Mathematically, one must solve two coupled partial dif- 
ferential equations 



d t q = d x [D(q) d x q] - d x [a{q) d x p] 



(22) 



and 

dtp = -D(q)d xxP - ~a'(q)(d x pf , (23) 

for the density field q(x, t) and the conjugate momentum 
field p(x, t). Here and in the following the prime denotes 
derivative with respect to the argument. Solutions with 
p(x, t) = are called relaxation solutions. For the relax- 
ation solutions Eq. (|23|) is satisfied, and Eq. (|22|) reduces 
to the hydrodynamic equation ([2]), so q{x,t) = p(x,t). 
Solutions with p(x, t) ^ are called activation solutions, 
here q(x, t) ^ p(x, t). Equations (|2"2"]l and (|23p are Hamil- 
tonian, 



d t q = SH/5p . 
with the Hamiltonian 



dtp = -SH/5q . 



H[q(x,t),p(x,t)} 



where 



3i(q,p) = -D(q)d x qd x p- 



dx "K, 



1 2 

-a(q)(d x p) 



(24) 



(25) 



(26) 



For a given model, specified by D and a, Eqs. ([22]) and 
can describe large deviations of different quantities 
in different settings. The problem of statistics of current 
during time T, starting from a step-like density profile, 
is specified by certain boundary conditions in x and t. 
To begin with, by virtue of mass conservation, a given 
integrated current implies an integral constraint 



J 



dx [q(x,T) - q(x,0)]. 



(27) 



The boundary conditions in t are different for the 
quenched and annealed settings. (The term 'boundary' 
emphasizes here that these are conditions at the bound- 
aries of the time interval [0, T].) In the quenched setting, 
the initial condition for the density coincides with Eq. (JXJ) : 



q{x, t = 0) = P-0(-x) + p+0(x), 



(28) 



where 6(x) is the Heaviside step function. The conjugate 
momentum is constrained by the condition at t = T [3GJ 



p{x,t = T) = \0{x), 



(29) 



where the Lagrangian multiplier A = A(J) is fixed by 
Eq. (gTJ). Once q(xA) and p(x, t) are found for < t < T, 
one can calculate [30| the function p, that enters Eq. © 
for the long-time asymptote of the moment-generating 
function: 



Mqucnchcd 



= A 

1 

2 



dx[q(x,T) - q(x,0)] 



dt 



dxa(q)(d xP ) 2 . (30) 
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The first term in ^ que nched comes from the constraint 
(f2"T)) . whereas the second one is equal to —S(T), where 
S(T) is the mechanical action of the Hamiltonian sys- 
tem and (HI). Indeed, 

/•T poo 

S{T)= dt dx{pd t q-W). (31) 

JO J -oo 

Using Eqs. (|22[) and (f26|) and performing integration by 
parts in the spatial integral in pip , one can rewrite the 
action as 



S(T) 



1 



dt / dx a(q){d x pf 



(32) 



In the annealed setting, the boundary condition at the 
final time t — T is again given by Eq. (l29l) . The initial 
condition is now different from Eq. (|2"5)) , it involves both 
q and p [3(| : 

p{x,0) = \8(x) + 2 / dr— h^-. (33) 
ip(x,0) cr(r) 

Once q(x, t) and p(x, i) are found, the function p can be 
calculated [30] from the equation 

Mannealod = - 2 / cfe / rfr — — - [g(o!, 0) - T\ 

J-oo Jp(x.fi) <*\?) 



+ A / efe [g(x, 1) - q(x,0)] 
Jo 

I r T r°° 
- - dt dx a(q)(d x pf . 

JO J-oo 



(34) 



The second and third terms here are the same as in the 
quenched setting, except that q(x,t) and p(x,t) are dif- 
ferent. The first term is specific to the annealed setting: 
it describes the cost of creating the optimal initial con- 
dition for the given value of the current. 



III. VARIANCE IN THE QUENCHED CASE 

From now on, we will only consider a class of models 
where D — 1 (such as in all three examples in Table 1). 
Here the governing equations (|22"|) and (|23| become 



dtq = d xx q - d x [a(q)d x p] , 
d t p = -d xx p- W (q)(d x pf 



(35a) 
(35b) 



The exact solution of Eqs. (|35a[) and (|35b[) . subject to 
the boundary conditions (|28|) and (|29|). is unknown ex- 
cept for the non-interacting random walkers, when o~(q) 
is proportional to q. Fortunately, the variance of current 
can be found for arbitrary a(q) by using a perturbation 
expansion over a hydrodynamic solution (qo,Po) = {p, 0), 
where 

( *\ P- + P+ , P+ - P- J X \ f OR\ 

P (x, t) = + erf f \ (36) 



solves Eq. <j2j) , with D = 1, subject to the initial condition 
(|28[). The Lagrangian multiplier A plays the role of a 
small parameter in this expansion. As one can justify 
a posteriori, a small A implies a small deviation of the 
current from its average value (J). We expand 



q = qo + Xqi + A q 2 + 
p = Xpi + X 2 p2 + 



(37a) 
(37b) 



and plug these expansions into Eqs. (I35aj) and (|35b|) . In 
the zeroth order we recover (qo,Po) — (/?, 0). The first- 
order equations are 



(d t - d xx )qi 
dtPi 



-d x [a(p)d x pi], 
-d xx pi. 



(38a) 
(38b) 



The boundary conditions for qi and p\ follow from 
Eqs. (HU and (|29]l: 

(39) 



g x (a:,i = Q)=Q ! Pi(x, t = T) = 9{x). 



Solving the anti-diffusion equation (j38bl) with the bound- 
ary condition (l39l) for p\ , we obtain 



(40) 



. . 1 1 
Pi{x,t) =2 + 2 



Now we need to solve Eq. (|38a|) : a diffusion equation with 
a known source term. The form of equation suggests to 
seek q\ as a gradient: 



qi = -d x tp. 
The potential ip satisfies the equation 

(dt - d xx )ip = F, 
where F — a(p)d x p\, and 
1 



d x p\ 



VMT - 1) 

The solution of Eq. dHJ is 



exp 



4(T - t) 



ip(x,t) 



In particular, 



i>(Q,T) 



dr 



F(y,r) 
dy — z exp 

, y/Mt-T) 



[x-v) 



(41) 
(42) 

(43) 

2 1 



4(t-T) 



dt 



dt 



°° , <r(p)d xPl 

ax — = exp 

-oo 



3C- 



VMT - 1) 

dxa(p)(d x pi) 2 . 



4(T - 1) 
(44) 



The function // from Eq. (30) becomes 

A 2 

Mquenched = A( J) + A 2 7^(0, T) - — ^(0, T) + 



A(J) 



-V(0,T) + O(A 3 ). 



(45) 
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Using ^ = A(J) + iA 2 (J 2 ) c + . 
the variance: 



(J 2 ) C = 0(O,T) 



dt 



[see Eq. (|7J|], we extract 



dxa(p)(d xPl ) 2 . (46) 



This result holds, for the quenched setting, for all diffu- 
sion processes with D(p) — 1 and arbitrary <r{p). 

The T-dependence of ( J 2 ) c can be easily extracted via 
the transformation t — > t/T and x — » x/s/T which re- 
duces Eq. (g6) to 



r ) n = vvt 



V = dt 
Jo 



dxa{p)(d xPl ) 2 . (47) 



Here p(x, t) is still given by (j3"6")l . whereas d x p\ is obtained 
from (|4"3"| by setting T = 1. Plugging these expressions 
into (j4"?| yields the announced result <j5j> - 

The variance (J 2 ) c corresponds to a Gaussian asymp- 
totic of the current distribution: 



P(J,T) 



1 



: exp 



%P)c 



(48) 



One can see from Eq. (|47p that the variance (J 2 ) c de- 
pends on the model only through cr(p). In simple models 
a(p) is a low-degree polynomial, see Table I. Consider 
now a more general case when a(p) admits a representa- 
tion 



n>0 



(The zeroth term in the series (j49|) vanishes, Aq 
since tr(0) = 0.) Combining (fMf and (|4*^|) we get 



quenched 



- E 

4-7T 



4?r ^ 2™ Uy p 



(49) 
- 0, 

(50) 



where d = p + — p_ , s = p+ + p_ , and 



t 



dx exp 



a;' 
2t 



erf 



The spatial integrals in E p vanish when p is odd, while for 
even p one finds E — \/8tt ,E2 = (3 — 2y2) \/87r, etc. 
The knowledge of E„ with p < 3 suffices to determine the 
variance for the 3-parameter class of models 

a = A lP + A 2 p 2 + A 3 p 3 . 

Here one obtains 

f 



^quenched 



8V27T 

3-2V2 

8V^ 



(4Ai s + 2A 2 s 2 + A 3 s 3 ) 
(2A 2 d 2 + 3A 3 d 2 s) . 



(51) 



We have also calculated the variance in the quenched 
setting from fluctuating hydrodynamics, see Appendix B. 



IV. VARIANCE IN THE ANNEALED CASE 

In the annealed case, the calculations are very similar, 
albeit somewhat more cumbersome. Employing the same 
perturbation expansion (|37a[) and (|37b[) . we recast the 
initial condition (1331) into 



p(x, 0) = A 
which yields 



(*) + 



Pl (x,o) = e(x) 



2giM) 
a[p(x,0)} 

291(1,0) 



0(A 2 ), 



(52) 



a[p(x,0)Y 

On the other hand, Eq. (|4U|) is still valid, and therefore 



Pl (x,0) = ± + ±E(x), £(x)=erf(| 



(53) 



(The T-dependence here is the same as in the quenched 
case, so we set T = 1.) Comparing Eqs. ({52} and (|53|) . 
we can deduce the initial condition on qy. 

^=\AT\fuVu x< >o < m) 

4 [<r\P+) L^W - 1] j z > 0. 

As in the quenched case, we employ a gradient represen- 
tation (|4"Tj) and find that the potential -0 satisfies the same 
inhomogeneous diffusion equation (|42[) . In the quenched 
case we had qi(x,0) = that led to ip(x, 0) = 0. In the 
annealed case the initial condition (|54|) leads to a non- 
trivial initial condition for ip: 



V>(x,0) = 0(-x)tr(p_) 



+ 0(a:)<7(^ + ) 



1 



£(*) 



1 



20F 

e- 2 /4 



4 

e(a:) - 1 



20? 



(55) 



where we have demanded that ip(x,0) be continuous 
at x — and chosen the arbitrary constant so that 
-0(0, 0) = 0. We now plug the A-expansions into Eq. ([34]) 
for ^annealed and obtain 



Manncalcd 



[qi{x,o)] 



dx ■ 

2„ 



A(J> +A 2 0(O,1) 
^2 ,1 

2 



dt 

•/ - X 



dxa(p)(d x pi) 



(56) 



with the same average current (J) as in the quenched 
setting. The variance is again extracted by using the 
A(J> + ±A 2 (J 2 ) C + . . .. The result is 



expansion p 

^annealed — 



dxa(p^)[£(x) + lf 



dx (j{p+) [£(%) — l] - " 



20(0,1) 



dt 



dxa(p)(d xP iY 
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After some algebra we find 



SSEP 



/>1 pOO 

^(0,1) = dt dxa{p){d xPl ) 2 

JO J-oo 

where 0) is given by (|55p . Combining everything and 
evaluating integrals, we arrive at the announced result ([6]) 
for the variance. 



V. EXAMPLES 

We now specialize the results to the three well-known 
models presented in Table f. Prior to that, however, 
we consider, for an arbitrary er(p), the symmetric case 
9- = P+ = P- 



A. Symmetric case 

For p- = p + = p the expression a(p) can be taken out 
of the integral in Eqs. ([5]) and ([6]), and we arrive at 



^quenched 
^annealed 



'2tt 



(57a) 
(57b) 



Hence in the symmetric case Kmnealed = >/2 Trenched 
independently of p and for arbitrary a(p). In the par- 
ticular case of SSEP, Eq. (|57bp has been known for a 
long time, see Ref. |38| and references therein, whereas 
Eq. (I57a[) (a gain, for the SSEP) has been obtained only 
recently j38j]. Derrida and Gerschenfeld [3(| noticed that, 
in models obeying the particle-hole symmetry, there is a 
symmetry relation between the optimal profiles in the 
quenched and annealed cases. For the SSEP, this rela- 
tion leads to 

/^annealed 

(A, 1/2) = v/2 /^■quenched 

(A, 1/2), (58) 

for /?_ = p + = 1/2 [3(J. We emphasize that relation 

Vanne aled = V% Vqucnchcd, which follows from EqS. (|57aj) 

and (|57bp . is valid for any p and it does not require the 
particle-hole symmetry, although we have derived it for 
only the variances. 



B. RWs 

For random walkers the variance is linear in the densi- 
ties p- and p+: 

tt P+ + P— Kinncaled rx Ircw 

('quenched = J== — , T7 = VA \pM) 

V27T V quenched 



Specializing Eq. (J5TJ) to A\ =2,A 2 = -2, and A 3 = 0, 
we find the variance for the SSEP in the quenched setting: 



2irV c 



quenched 



p+ + p- 
3-2V2 



{P+ + P-Y 



(p+-p-Y 



(60) 



To our knowledge, this result is new. In the annealed 
setting we get 



Vanncalcd = 9- + P+ ~ P- ~ P+ 



(61) 



in agreement with Eq. (|13[) . In the extreme anti- 
symmetric case (/?_,/?_!_) = (p, 0) the variance reads 



^quenched 
^annealed 



V2 



(62) 



One can see that Vquenchcd < Knncaicd for all < p < 1; 
the equality occurs only when p = or p = 1 . 



D. KMP 

For the KMP model we have A 2 = 4 and A\ = A 3 = 0. 
As a result, 



quenched 
^annealed 



( P+ +p_) 2 + (3-2V2)( P+ -p_) s 
±p-p+ 2(p--p+) 2 



'2tt 



The expression for V annoa i c d coincides with C 2 from 
Eq. (fT6|). the expression for lq Ue nched is new. 



VI. LARGE DEVIATIONS: A NUMERICAL 
SOLUTION 

What happens when A is not small or, in other words, 
the current J(T) is not close to the average current (J)? 
As we cannot solve the MFT equations (|22|) and (|23[) an- 
alytically, we resort to a numerical solution. Bunin et 



in agreement with Ref. [3C 



al. [32j have recently developed a numerical algorithm 
based on minimization of the mechanical action. In their 
algorithm, the boundary conditions in time involve the 
knowledge of the density profiles at some initial and final 
times. In the context of integrated current fluctuations, 
one needs an algorithm that would deal with a boundary 
condition on the momentum at t — T. Fortunately, clas- 
sical field-theoretic Hamilton equations have previously 
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appeared in many different contexts. It is hardly sur- 
prising, therefore, that an efficient and simple numerical 
algorithm, with the required type of boundary condition 
at t = T, already exists. It was originally suggested by 
Chcrnykh and Stepanov 31] for evaluating the probabil- 
ity distribution of large negative velocity gradients in the 
Burgers turbulence. Later on it was employed by Elgart 
and Kamenev [ljj and Meerson and Sasorov [l8| for eval- 
uating the mean time to extinction in finite-size lattice 
gas systems involving random walk and on-site reactions. 

The algorithm iterates the diffusion-type Eq. (22) for- 
ward in time and the anti-diffusion-type Eq. (|23|) back- 
ward in time. Consider first the quenched case. In a 
simple version of the algorithm, each iteration of q(x, t) 
starts at t = from the initial condition (|28[) and solves 
Eq. (|22l) forward in time until time t = 1 is reached. In 
this calculation the previous iteration for p(x, t) is used. 
Then Eq. (|23|) for p is solved backward in time starting, 
at t = 1, from p(x, 1) = X6(x) [see Eq. (|2"5|) ]. and contin- 
uing until t — 0. Here the previous iteration for q(x, t) is 
used. The very first iteration for p is simply the desired 
final state p(x, 1) = X6(x). 

Unfortunately, the simple version of the algorithm suf- 
fers from a numerical instability: after an initial tran- 
sient, the numerical solution alternates between two dif- 
ferent sets of q and p, instead of con verg ing to a unique 
(q,p) solution [39j]. Similarly to Ref. [3l| we suppressed 
this instability by replacing p, in the iterations for q, by 
a linear combination of the values of p obtained in two 
previous iterations. Similarly, we replaced q, in the itera- 
tions for p, by the same linear combination of the values 
of q obtained in two previous iterations. The relative 
weights of these two values of p and q (the coefficients of 
the linear combination) must sum up to unity. In the ex- 
amples shown below we chose the previous iteration with 
weight 0.75 and the iteration before previous with weight 
0.25. The first two iterations of q and p are performed 
with the simple version of the algorithm. 

For the annealed setting we have to use the initial con- 
dition (|33p which involves both q and p which are a priori 
unknown. Therefore, in the very first iteration we solve 
Eq. (f2"2"|) for q forward in time, starting from a "wrong" 
(quenched) initial condition, Eq. (|28l) . Then, after solv- 
ing Eq. ([23]) for p backward in time until t = 0, we de- 
termine q(x, 0) from Eq. (f3"3"]) . feed it into the forward- in- 
time solution for q, and continue iterations. Here too we 
used, starting from the third iteration, the values from 
two previous iterations to suppress the numerical insta- 
bility. 

We implemented this algorithm in Mathematica. We 
worked with a finite-size system, \x\ < L/2 and imposed 
the boundary conditions q{—L/2,t) = p_, q(L/2,t) — 
p+, p(-L/2,t) = 0, and p(L/2,t) = X. The step- 
functions entering the boundary conditions at t = and 
t = 1 were smoothed a bit. The iterations converge 
very rapidly. Having computed q(x,t) and p(x,t), one 
can evaluate the large-deviation function p by numeri- 
cally evaluating the integrals in Eqs. (|30p and (|M]l for 




FIG. 2: (Color online) A numerically computed optimal path 
for the SSEP with a step-like initial density profile (p_ = 0.6, 
p+ = 0.2) in the annealed setting. Shown are the particle 
density q(x, t) (upper panel) and the conjugate momentum 
density p[x,t) (lower panel) at rescaled time moments t = 
(solid line), 1/3 (dashed line), 2/3 (dotted line) and 1 (dash- 
dotted line). The Lagrangian multiplier A = 4 corresponds to 
the rescaled current J = p'(X)\\ = 4 = 1.118..., where p(A) 
is taken from Eq. (|10|) . For comparison, the rescaled average 
current is (J) = (p_ — p + )/ v / 7r = 0.2256 .... 



the quenched and annealed settings, respectively. 

Figure [5] shows an example of numerically found op- 
timal path q(x,t), and the corresponding p(x, t), for the 
SSEP in the annealed setting. In this case the function 
p is known, see Eq. ([101 . A = 4 corresponds to a positive 
current about five times greater than the average current. 
The optimal initial profile is such as to facilitate hydrody- 
namic transport that does not cost action. The optimal 
fluctuation grows towards t = 1 , when the hydrodynamic 
flow weakens. In this example our numerically computed 
values of the function p and rescaled integrated current 
[using Eqs. (f34|) and ([3]), respectively] are within 3 and 
4 per cent, respectively, of their theoretical values. The 
numerically found values of p for A = 4 and —4 are shown 
in Fig. [1] alongside with the exact and approximate ana- 
lytical results for p(X). 

Figures |3] to [3] show three examples of numerically 
found optimal paths q(x, t), p(x,t) for the SSEP in the 
quenched setting. Here no analytic results are available 
beyond the small fluctuations, see Eq. ([5]), except for the 
special cases of p_ = p + = 1/2 and p_ = 1, p + = 0. 
The case of A = 4 corresponds to a positive current a 
few times greater than the average current. Here too 
the optimal fluctuation grows toward t — 1 and facili- 
tates transport of the material from left to right. Notice 
a striking similarity between the p-profiles for A = 4 in 
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FIG. 3: (Color online) Same as in Fig. [2] but in the quenched 
setting. Here A = 4 corresponds to the rescaled current J ~ 
0.9, as found numerically. 



FIG. 4: (Color online) Same as in Fig. [3] but with A = — 4 
which corresponds to a negative rescaled current J ~ —0.3, 
as found numerically. 



the annealed and quenched settings, for which we do not 
have a good explanation. 

For A = —4 the optimal fluctuation reverses the cur- 
rent compared with the hydrodynamic flow. Finally, for 
A ~ —1.34 the integrated current is equal to zero. Here, 
after an initial release of material from left to right, the 
fluctuation pushes the material back. Still, as one can 
see from Fig. [5j the final density profile q(x, 1) is differ- 
ent from the initial profile. 

Finally, Fig. [6] shows the A-dependence of 
Mquenched/'V^ tnat we f° un d numerically in a moder- 
ate range of A for p_ = 0.6 and p + = 0.2. Also 
shown is the two-cumulant approximation fi = C±X + 

(1/2) VquenchedA 2 , with Ci from Eq. (JT2J and trenched 



from Eq. (f60j) . One can see that, as in the annealed set- 
ting, the two-cumulant approximation is quite accurate 
well beyond small A. That is, both in the annealed and 
in the quenched settings, deviations of the probability 
distribution P(J,T) from a Gaussian only occur in rela- 
tively far tails of the distribution. 



VII. CONCLUDING REMARKS 

We have investigated the long-time fluctuations of inte- 
grated current in diffusive lattice gases in one dimension, 
when the initial density is a step-like function. Our anal- 
ysis relies on the macroscopic fluctuation theory (MFT) 
13], more precisely on its implementation (30j which al- 
lows one to examine the fluctuations of the current in 
a non-stationary situation of a step-like initial density 
profile. For the quenched and annealed settings we have 




FIG. 5: (Color online) Same as in Figs. [3] and [4] but with 
A = —1.36 which corresponds to J ~ 0. The rescaled time 
moments are t = (solid line), 1/2 (dashed line), and 1 (dot- 
ted line). 



calculated the variance of the current fluctuations by de- 
veloping a perturbation theory around the noiseless hy- 
drodynamic solution. Our results for the variance hold 
for a whole family of lattice gas models which, at the 
coarse-grained level of the MFT, can be characterized 
by a constant diffusion coefficient and an arbitrary a(q). 
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A 

FIG. 6: (Color online). Mqu^Tchcd/v 7 ? versus A for p- = 0.6 
and p+ = 0.2. Circles: our numerical results. Dashed curve: 
the two-cumulant approximation CiA + (l/2) U) Ue nchedA 2 with 
Ci from Eq. (fT2"1) and V qucnc h c d from Eq. (pH) . 



sponding to an unusually large current, to a current that 
flows in the "wrong" direction, and to zero current. We 
have also computed numerically the large deviation func- 
tion fi(X) and observed that the two-cumulant approxi- 
mation (and, correspondingly, the Gaussian asymptotic 
of the integrated current distribution) remain quite accu- 
rate well beyond the small-A regime. These numerical re- 
sults are important because an analytical solution of the 
MFT equations, beyond small fluctuations, is presently 
unavailable except for non-interacting random walkers 

M- 

An important task for a future work is an asymptotic 
analysis of the MFT equations in the large current limits, 
\J\ 3> (J). Here the iteration algorithm may help in 
getting insight into the type of perturbation expansion 
needed for that. 



Particular examples of our results include the variance for 
the non-interacting random walkers, for the symmetric 
exclusion process and for the Kipnis-Marchioro-Presutti 
model. For the annealed setting these particular results 
agree with previous results. For the quenched setting 
these results are new. 

We have also investigated numerically the regime of 
large deviations of the current. Using the SSEP as a 
representative example, we have solved the MFT equa- 
tions by using the Chernykh-Stepanov numerical itera- 
tion algorithm. We have found the optimal paths corre- 
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Appendix A. Free energy from the MFT formalism 



relaxation trajectory. That is, the optimal fluctuation 
q(x,t) must obey the time-reversed version of Eq. @: 

d t q = -d x [D(q)d x q}. (A3) 

Combining this equation with Eq. ([22)1 . we obtain 

a(q)d x p = 2D(q)d x q + f(t), (A4) 

where f(t) = because of the boundary conditions in x. 
Now we introduce function F(q) that satisfies Eq. (|21|) . 
Integrating Eq. (|A4p over x yields 

p = F'(q)-F'(p), (A5) 

where we have demanded that p = at q — p. A straight- 
forward algebra shows that p — F'(q) — F'(p) also solves 
Eq. (|23| . This local relation between p and q implies com- 
plete integrability of the MFT problem for equilibrium. 
Furthermore, here "K = 0, and the action So becomes 

/0 />oo 
dt / dxpdtq 
-oo J — oo 

/oo />0 
dx dt[F'(q)-F\p)]d t q (A6) 
-oo J — oo 

which, upon integration over time, yields Eq. (| Al|) as 
expected. 

In non-equilibrium situations the local relation p = 
F'(q) — F'(p) breaks down, as it does not satisfy some or 
all of the boundary conditions. As a result, the MFT 
equations become, in general, non-integrable, whereas 
the mechanical action explicitly depends on the system 
dynamics at intermediate times. 



Here we show how Eq. (|2Tj) appears in the MFT for- 
malism. (See also Ref. [Tj] for a similar derivation in the 
particular case of the SSEP.) Consider a lattice gas in 
equilibrium at average density p = const. The probabil- 
ity of observing a given density profile q(x) is described 
by the equilibrium Boltzmann-Gibbs distribution: 

/>oo 

-ln7[q(x)]~ dx [F(q(x)) - F(p) - F'(p)(q(x) - p)] , 



FT! (A1) 
see e.g. Ref. [30(. In the MFT formalism, this expression 

should coincide with the mechanical action, 

/0 poo 
dt / dx(pd t q-J{), (A2) 
-oo J —oo 

calculated along the activation trajectory of the MFT 
equations and P3"|) obeying the following boundary 
conditions in time: q(x,t — —oo) — p and q(x,t = 0) = 
q(x) [HI- We will now see how it happens, and how F(q) 
emerges. 

It is very simple to find the activation trajectory here 
because, in equilibrium, it coincides with a time-reversed 



Appendix B. Quenched variance from fluctuating 
hydrodynamics 

Equation ([5]) for the variance in the quenched setting 
can be also obtained in the framework of fluctuating hy- 
drodynamics, by solving a linearized Langevin equation. 
Once D = 1 and o~{p) are known, the Langevin equa- 
tion for the fluctuating particle density field q(x,t) can 
be written as 0, 0] 



d t q = d xx q + d x \/(j(q)£(x, i) 



(Bl) 



Here ^{x, t) is a zero-average Gaussian noise, delta- 
correlated in space and in time: 

(ax,m(x 1 ,t 1 ))=6(x-x 1 )6(t-t 1 ), (B2) 

and the brackets denote ensemble averaging. Linearizing 
the Langevin equation (|B1|) around the hydrodynamic 
solution (|36|) , we write q = p + qi where qi <§; p. This 
yields 



d t qi - d xx qi = d x y/a(p)£(x, t) 



(B3) 
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Plugging qi = d x <fi, we can rewrite this equation as The variance of the current is, therefore, 



8^-8^= yMp)Z(x : t). (B4) 



The solution is (J }c = ( J dx J dy q 1 (x,T)q 1 (y,T) 



OO f DO 



f dt x T dy^^kMe~^{Bb) = ( r dx ^ dyd x cb(x,T)d y cj)(y,T)) 

Jo J-oo \/47r(t-ti) \Jo Jo I 

= (0 2 (O,T)). (B6) 



The current J(T) from Eq. (|27p can be represented as 



/■OO 

J(T) = (J) + / gi(af,T) dx. Combining Eqs. |B5} and (|B6|) . we obtain 

Jo 



„ 1 [ T f T f°° [°° x 2 V 2 . 

J )c = -r dt i dt i / dx / 4<T ~ 4l) ^ vVRm^R^)] t 2 )) 

47T Jo JO J-oo J-oo 



JO J —oo J —oo 

T r T r oo r oo 



^-J dhj ,11, I ,/.,' / <h,<->-T r :i , v ■a[ l , { :r.l l )] n [,, {!/ .l,) Sir- n)HU -I-,) 

^ pT j-OO — 2 (T_t) 



sy. "'L**— <*» < B7 > 



which coincides with Eq. (|46[) obtained, in the quenched 
setting, from the MFT. 



